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Abstract 
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I. INTRODUCTION 



In the early 1980s single-mode silica optical fibers were recognized as a privileged medium 
for experiments in quantum optics because they exhibit a well defined transverse mode and 
very low losses. The Kerr nonlinearity of silica has been used to produce non classical states 
of light. Squeezed state production in fibers was pioneered by Levenson et al. WW in 1985. 
Since then a large number of squeezing experiments have been performed with cw-light or 
short optical pulses (for a review see 3(). 

More recently, Fiorentino et al. demonstrated that optical fibers may also be used to 
produce twin-photons pairs 4]. This new kind of twin-photon source is well suited for fiber- 
optic quantum communication and quantum cryptography networks. In contrast with the 
traditional twin-photon sources, based on the down-conversion process, it avoids large 
coupling losses occurring when the pairs are launched into long distance communication 
fibers. The physical process used in to generate twin-photon pairs is a four-wave mix- 
ing (FWM) phase-matched by the interplay of the optical Kerr effect and the chromatic 
dispersion, called (scalar) modulation instability (MI). This process can be basically under- 
stood as the destruction of two pump photons at frequency u and the creation of Stokes 
(red-shifted) and anti-Stokes (blue-shifted) photons at frequencies u> s and u a satisfying the 
energy conservation relation 2u; = ^ s + cu a . In the time domain, the beating of the pump, 
signal and idler waves produces a fast modulation of the pump envelope. 

MI is a spontaneous phenomena which can be described in the framework of classical 
nonlinear optics In this classical framework, one usually considers that MI is induced 
by some incoherent noise initially present on the pump wave. Twin-photon pair production 
requires however a very low input noise power to be efficient (otherwise the twin-photons 
are buried in the background photon noise). Such a regime is dominated by vacuum fluctu- 
ations and requires a quantized field theory to be described properly. A small-perturbation 
approach of this problem (non depleted pump approximation) shows that vacuum fluc- 
tuations can induce the MI even in the absence of any classical input noise. This is an ideal 
situation for twin-photon pairs production, but it never occurs in a real life experiment. In 
practice, classical input noise and vacuum fluctuations compete for inducing MI. 

Understanding the dynamical development of MI from a quantum point of view is an 
important issue for the design of fiber-optics twin-photon pairs sources. We present an 
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unified approach to this problem based on the stochastic nonlinear Schrddinger equations 
(SNLSE) pi 0, 0, [ill [ill [l^. This formalism is equivalent to quantum-field Heisenberg 
equations but has two main advantages. First, it is very suitable for numerical simulations 
of complex situations where classical noise and vacuum fluctuations act together. Second the 
correspondence principle of quantum physics, i.e. the transition from quantum to classical 
world, looks very natural in the SNLSE formalism. 

The MI observed in standard non-birefringent single-mode silica fibers is often referred 
to as scalar modulation instability (S-MI) because polarization of light plays no role in this 
process. In birefringent fibers different kinds of MI may appear because of the interplay 
of nonlinearity, chromatic dispersion and birefringence. This case will be referred as vector 
modulation instability (V-MI). Twin-photon pairs sources based on V-MI have not been 
demonstrated yet, in contrast with S-MI based sources j^J. In this article, we will however 
address both issues because we recognize that V-MI twin-photon pairs sources could be 
more practicable than S-MI based ones. The theoretical and computational tools that we 
developed apply to all kind of MI, including the low-birefringence, high-birefringence and 
scalar limits. 

This article is divided into four sections, beginning with the present introduction. The 
second section consists of a review of the approach based on a perturbation analysis around 
the steady state solution. This is the simplest method of approaching the problem of MI. 
It allows us to set the stage for the more sophisticated approach we then develop and to 
make contact with earlier work in the field. Then in Sec. IIHI we exhibit the SNLSE for 
scalar and vector MI. Contrary to the approach of Sec. |TI] these equations are not based 
on a perturbation analysis and are valid in the strongly non linear regime. They describe 
both MI originating from classical noise and from vacuum fluctuations (quantum noise). 
The SNLSE we obtain generalise the earlier resnlts of fl fl Q R Q. In particnlar the 



SNLSE we obtain does not require the birefringence to be either small or large as in earlier 

nn 

work PJJ, llil, but are valid for all values of the birefringence. For the interested reader 
we present a self contained derivation of these equations in the Appendix. In Sec. IIVI we 
use the split-step Fourier method to integrate the SNLSE derived in Sec. IIHI and illustrate 
our algorithm on the cases of scalar MI and of vector MI both for weak, intermediate, 
and strong birefringence. In order to interpret the results of the numerical integration it is 
essential to introduce the notion of mode. This allows us to compare the numerical results 
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and the analytical solutions derived from perturbation theory. We also compare in detail 
the characteristics of MI induced by classical noise and by quantum noise. In a companion 
article we shall show that our numerical results are in very good quantitative agreement 
with experimental results. 

II. SCALAR MODULATION INSTABILITY: PERTURBATION ANALYSIS 

As pointed out in the introduction several kind of MI can occur in optical fibers. From 
a general point of view, MI occurs when the continuous wave (steady-state) propagation 
is unstable. The most straightforward way to get some insight on how the MI develops 
is to perform a small-perturbation analysis of the steady-state. However this method has 
limitations: it cannot address the strongly non linear regime and simple analytical formulas 
cannot be always obtained by this method. 

In this section we illustrate the perturbation method in the case of S-MI, which is the 
simplest. S-MI occurs in isotropic fibers in the anomalous dispersion regime. In what follows 
we will review the most important aspects of this approach, focussing on the comparison 
between the classical and quantum descriptions, and on the limitations of this approach. 
In Sec. IIVI we will compare the analytical formulas for S-MI derived from the perturbation 
analysis to numerical results from the SNLSE derived in Sec. 11111 The discussion of MI in 
birefringent fibers is also postponed to Sec IIHI and Sec. IIV1 

A. Classical description 

In an isotropic single-mode fiber, the electric field may be written 

E(r, t) = F(x, y)A(z, t) exp[if3 z — u t]Si + c.c. (1) 

where x is a unit vector orthogonal to the fiber axis (z-axis), ujq the carrier angular frequency 
and (3q = (3(uJo) the associated propagation wave number (modal propagation constant). 
F(x,y) stands for the mode profile function and A(z, t) for the field envelope. The field is 
supposed to be polarized linearly. This is not restrictive because of the isotropy assumption. 
The complex envelope A evolves according to the nonlinear Shrddinger equation (NLSE) that 
can be established from Maxwell's electromagnetic theory [5|. If the envelope is normalized 
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in such way that |A(z,i)| 2 is equal to the instantaneous power flowing through the plane 
z = constant at time t, the NLSE is obtained in the following form: 

OA 1 dA .fad 2 A , . A 

Here v g = (dfi/duj)* 1 is the group velocity of the wave, fa = d 2 f3/du 2 the group- velocity 
dispersion (GVD) parameter and 7 is fiber nonlinearity parameter defined as 



4e ngc 2 Aeff 



7 = T. ^2^2 A ( 3 ) 



where n$ is the fiber mean linear index of refraction, A e s the mode effective area and Xxxxx 
is the diagonal element of the fiber x nonlinearity tensor. 

S-MI is observed when a continuous wave (cw) or a quasi-cw optical pulse is launched 
into the fiber. In the cw-case, when an optical power Pq at frequency ujq i n injected, a 
first order perturbation analysis of Eq. (j2J) shows that the cw steady-state solution A st (z) = 
y^exp (i'-fPoz) becomes instable in the anomalous dispersion regime (fa < 0). The insta- 
bility manifests itself by a parametric gain at frequencies ujq ± VI with < < 2W r yPo/\fa\. 
The maximal gain g = 27 P L occurs for f2 = ■\/2 r yP /\fa\ = f2 max 5]. Noise at these fre- 
quencies is strongly amplified. As a result, the optical spectrum at the fiber output exhibits 
two sidebands at frequencies u ± f2 max - This analysis supposes that the pump power re- 
mains constant (undepleted pump approximation) and a cw regime. However, numerical 
simulations of Eq. (J2J) show that the above formula also hold for quasi-cw (for instance 
nanosecond) optical pulses (see Sec. II V|) . 

Sideband generation can be interpreted as a FWM process phase-matched by the interplay 
of dispersion and Kerr nonlinearity. From this point of view, the MI process can be seen 
as the destruction of two pump photons at frequency ujq followed by the creation of Stokes 
and anti-Stokes photons at frequencies u> a = uo — fi max and u a = u + fi max respectively. 
In the cw-case and non depleted pump approximation, the output power spectral density 
at frequencies u s and u a can be computed using analytical formulas |3| for given initial 
conditions. These formulas also hold for quasi-cw pulses provided one replace Stokes and 
anti-Stokes power spectral densities by the number n s and n a of Stokes and anti-Stokes 
photons located in the same temporal mode as the pump wave (see Sec. HVBI for a detailed 
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analysis of this issue). For example, assuming an incoherent initial noise, one finds 

n s (L) = n s {0) cosh 2 ( 7 P L) + n a (0) sinh 2 ( 7 P L), (4a) 
n a {L) = n a (0) cosh 2 ( 7 P L) + n a (0) sinh 2 ( 7 P L), (4b) 

where L is the propagation distance. These equations show that classical incoherent noise 
is amplified exponentially with gain g = 2 7 P L. 

Noise induced S-MI was first demonstrated experimentally by Tai et al. in 1986 Q|. 
This experiment confirmed two main predictions of the classical theory: the amplification of 
Stokes and anti-Stokes photons and the power-dependence of their frequency shift. However 
the classical theory of S-MI only holds when the initial number of noise photons is high or in a 
stimulated regime when a coherent probe pulse at Stokes or anti-Stokes frequency is injected 
together with the pump pulse Q|. Eqs. (J3J are unable to explain ab nihilo generation of 



Stokes and anti-Stokes twin-photons reported in 4J. Neither are they valid when the mean 
photon numbers n s (0) and n a (0) are of the order of one or lower. In this regime vacuum 
fluctuations play a central role and field quantization is required. 



B. Quantum description 

To take into account vacuum fluctuations, fields must be quantized. The quantum coun- 
erpart of the NLSE (J2J) is known as the quantum nonlinear Shrodinger equation (QNLSE) 
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where A is the quantum operator corresponding to the field envelope and A^ its hermitic 
conjugated. A and A* satisfy the bosonic equal-space commutation rule: 

A(z, t), ti{z, t')] = hw 5(t - £')■ (6) 

The normalization constant hwo stems from the normalization chosen for the field envelope. 
In the quantum propagation theory, the expectation value of the optical power flowing 
through the plane z = constant at time t is given by (A'(z,t)A(z,t)). In a one dimensional 
system, space and time play a symmetrical role. The dynamics can be described either 
in terms of spatial wave-packet evolution (evolution picture) or in terms of temporal wave- 
packet propagation (propagation picture). The first picture is the must common in quantum 



field theory and leads to equal-time commutation rules between the envelope fields A and 
A 1 . In this section however, we chose to work in the propagation picture (which is usual in 
nonlinear optics) in order to get a closer correspondence between quantum (jSJ) and classical 
^ equations. In this case equal-space commutation rules © must be imposed to A and 
A* |2l|. (In contrast, the evolution picture is used in the Appendix D to derive the stochastic 
equations of Sec. 11111 ) 

In the cw regime, the basics of the quantum theory of S-MI can be understood by per- 
forming a first-order perturbation analysis of the steady state solution of Eq. (J3J). Be- 
cause the pump field contains a large number of photons one can treat it as a classical 
coherent wave. Using this approximation the steady-state solution is just the classical one: 
A s t(z) = -\/Pd ex P (ijPoz)- The disturbed field can be written as: 

A(z,t) = A st (z) + a(z,t). (7) 

The disturbance operator a(z, t) is defined by Eq. (J7J). It satisfies the same commutation rule 
as A(z, t): [a(z, t), aJ(z, t')\ = hu 5(t — t'). Injecting the ansatz ((7j) into Eq. (jSJ) one obtains 
a propagation equation for the disturbance a(z,t). Supposing that the disturbance is small 
(non depleted pump approximation) this equation can be linearized and solved analytically 
in the Fourier domain jfjj]. 

Using this method one finds that the quantum and classical theory predict the same 
frequency dependence of the parametric gain. In contrast the quantum and classical theory 
differ in that they do not predict the same growth law for Stokes and anti-Stokes photon 
numbers. One easily finds that the quantum counterparts of Eqs. (J3J are 

n s (L) = n s (0)cosh 2 ( 7 P L) 

+ K(0) + l]sinh 2 ( 7 P L), (8a) 
n a {L) = n a (0)cosh 2 ( 7 P L) 

+ K(0) + l]sinh 2 ( 7 P L), (8b) 

where n s and n a stand now for the expectation values of the photon number operators: 
Ui = {fij), i = s,a. These equations show that S-MI can be observed even in absence of any 
classical input noise. Setting n s (0) = n a (0) = in Eqs. (jHJ) gives the number of Stokes and 
anti-Stokes photons produced by the sole action of vacuum fluctuations: 

n s (L) = n a (L) = sinh 2 ( 7 P L). (9) 
7 



Eqs. (JHJ) hold only for perfectly phase-matched photons at frequencies uj s and u a and uncor- 
related initial noise but can be easily generalized to get round these restrictions. 

Although Eqs. (JSJ) and (JHJ) give a good insight into the physics of vacuum-fluctuations 
induced S-MI and photon pairs generation, they are not suitable for quantitative predictions 
in the pulsed regime. This is because the effective value of the pump power P depends on 
pump pulse shape, duration, and spectral width. Furthermore the energy spectral density of 
Stokes and anti-Stokes waves deduced from is highly dependent on the precise definition 
of modes. This will be discussed in Sec. IIVBI In the next sections we will show how to get 
around these difficulties by introducing the SNLSE and solving it numerically. 



III. SCALAR AND VECTOR MODULATION 
INSTABILITIES: STOCHASTIC EQUATIONS 

One can go beyond the perturbation analysis of S-MI by solving the QNLSE (JSJ) numer- 
ically. Such a plan could seem cumbersome because Eq. (JSJ) is a field-operator equation. 
The problem can be bypassed by converting the operator equation (JSJ) into c-number equa- 
tions. This can be performed by choosing a representation for the electromagnetic field. In 
this article, we will use the positive P- representation (P^) introduced by Drummond and 
Gardiner |22|. The c-number equations obtained in this way are not standard deterministic 
partial derivative equations but stochastic (Langevin-type) ones. 

Using the P^ representation, it can be shown 0, 0, E| that the QNLSE (JSJ) is equivalent 
to the following set of stochastic equations: 



dA 1 OA 3 d A 

~dz~ + ~~dt = - J y^ + ^[ it (^) i (^)] A (^) + &CM) A(z,t), (10a) 

~dz~ + ~~df = +l 2~dP _Z7[ ( Z ^) A ( Z ^ A + V-htwo Uz,t) A\z,t). (10b) 

Eqs. (1 10b ) and ()10b) look like the classical NLSE (J2J) and its complex conjugated, except 
for the last terms which accounts for vacuum fluctuations. These last terms contains two 
independent zero-mean Gaussian white noise random fields (i(z, t) and £2(2, t) characterized 
by the following second order moments: 

(C k (z,t)d(z',t')) = 5 kl 5(z- z')5(t-t'), (11) 
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with (k,l) G {1,2} 2 . Because the random fields (i and (2 are not complex conjugated of 
each other, the envelope fields A and are only complex conjugated "in mean" and have 
to be treated as different mathematical objects. 

Eqs. (|1U|) can be solved on a computer. The numerical methods used for this task will be 
briefly explained in Sec. II V Al Note that solving Eqs. (|1(J|) gives a single realization of the 
stochastic process. In order to calculate the expectation value of a quantum observable, a 
statistical average on many realizations is required. Thus one generates a large number of 
realizations (A[ n ](z, t), At,(z, t)), n = 1, . . . , N. In order to calculate the expectation value 
of a quantum observable one then carries out a statistical average over the many realizations. 
For example, the expectation value of the energy spectral density of a pulse is given by: 



where X(z, ft) = J™ X(z,t) e lQ,t dt designates the Fourier transform of the field X(z,t) 
and Q is the detuning from the pump angular frequency ujq. 

In practice, a few realizations are enough when the number of photons per mode is 
high. However in a regime dominated by vacuum fluctuations hundreds of realizations are 
typically required to get precise values. A comparison between Stokes and anti-Stokes photon 
production predicted by the SNLSE ()10|) and the analytical formulas (jSJ) will be presented 
in sections ITVBl and WU 

The V-MI occurs in birefringent fibers. These fibers are characterized by different propa- 
gation constants (3 Qx and (3 0y and different group velocities v gx and v gy for the x- and y-axis 
polarized modes. The numerical study of vacuum fluctuations induced V-MI requires an ex- 
tension of Eqs. (fTUj) that takes into account the phase mismatch parameter A/3 = (3q x — (3q v 
as well as the group- velocity mismatch A/3i = l/v gx — l/v gy . Such extensions have been 
established in earlier works for the low-birefringence Q| and the hight-birefringence Q| 
limits. We generalized these results for an arbitrary level of birefringence and obtained the 
following stochastic coupled nonlinear Schrodinger equations: 




(12) 
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where (Aj;,4) an d (4' 4) are stochastic envelope fields associated to the x- and y-axis 
modes respectively, and B = Xxyyx/Xxxxx is a parameter that measures the strength of the 
nonlinear coupling between the x and y components. Its value lies between and 1 and 
depends on the nonlinearity mechanism. For silica fiber, we can set 5 = 1/3 because the Kerr 
non linearity has principally an electronic origin. Four independent gaussian random fields 
(k(z,t) are needed to reproduce the effect of vacuum fluctuations. They are characterized 
by the second order moments (|lljl. as in the scalar case, with (k,l) G {1,2,3,4} 2 . The 
demonstration of this set of equations is outlined in the Appendix B 

Note that, in contrast to the scalar perturbation analysis does not lead to sim- 

ple analytical formulas for vacuum-fluctuations induced V-MI. For this reason Eqs. ([13^1 
constitute a valuable theoretical tool. 



IV. NUMERICAL INTEGRATION OF THE SNLSE 



A. Sample Spectra 



We have developed a method for integrating numerically the SNLSE in the case where 



the pump is representee . 
Fourier (SSF) method 



?y a pulse of finite duration. Our method is based on the split-step 
5[. We have had to generalize the method in two ways. First of all 



the stochastic noise is modeled by including a noise term at each propagation step. Second 
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in the case of the V-MI we must alternate not only between the time and Fourier domain, 
but also between the linear and circular polarization basis. Switching from time to Fourier 
domain is the basics of the SSF method: it permits to handle the time-derivatives in a 
simple way. Similarly, it turns out that, whereas the terms with time- derivatives are easier 
to handle in the linear polarization basis, the 7-dependent terms are better managed in the 
circular polarization basis. 

The quantum noise in the SNLSE (JT5Jl contains four independent real zero-mean Gaussian 
white noise functions Ck( z ) t) characterized by Eq. pip. In the numerical method, time and 
space are discretized with respective discretization steps r and h. So each family of noise 
functions Ck{^,t) becomes a finite number of random variables Ck chosen according to 
a zero- mean Gaussian law of variance The variance value is imposed by the normaliza- 
tion condition pip . The matrixes Cfc^z^t] (k = 1, 4) define the stochastic path of each 
realization. In contrast, when we will study the effect of classical noise we will add it once, 
at the beginning of the pulse propagation, to the spectral distribution of the signal A(0, fl). 

As we have previously indicated only the expectation values of observables have a phys- 
ical meaning in the stochastic equations. From a numerical point of view this means that 
one must average the calculated quantities over several realizations of the stochastic path, 
and/or the classical input noise. Usually, averaging over a hundred of realizations gives an 
uncertainty on the numerical results less than 1 dB in the non-zero gain frequency-range. Fi- 
nally we note that including the stochastic terms do not increase significantly the numerical 
complexity of a single realization. 

Some sample spectra obtained using our algorithm are presented in Fig. ^ in the cases 
of S-MI (Fig. [T^,), low birefringence V-MI (Fig. HJ)), and high birefringence V-MI (Fig. [TJ;). 
The energy spectral density is plotted versus the frequency detuning from the pump. The 
physical parameters used in these simulations are listed in Table |l] In every simulation, the 
pump wave has been supposed to be an unchirped linearly polarized Gaussian pulse with 
peak power P Q and full- width-half-maximum duration Tfwhm- No classical noise has been 
added, and an average over 50 realizations of the stochastic process has been performed. 
The frequency detuning of the side bands agrees with linear perturbation theory. (In order 
to make the comparison easier we have indicated in Fig. ^i-c the angular frequency shifts 
^max at which maximum gain is expected on the basis of the linear perturbation theory.) 

Moreover the SNLSE predict quantitatively the effect of vacuum fluctuations on the 
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TABLE I: Simulation parameters for Fig. ^i-c. 

evolution of the energy spectral density of the electromagnetic field. In Sec. IIVBI we will 
show that this evolution is also in very good agreement with linear perturbation theory when 
the number of Stokes and anti-Stokes photons generated is small enough. When the side 
bands are well-developed, the perturbation theory fails to predict correct values of Se- In 
contrast, the SNLSE algorithm still gives accurate results. In this limit numerical results 
can be easily confronted to experimental data. In Q], we report an experiment on high 
birefringence spontaneous V-MI in the anomalous dispersion regime that shows that the 
theoretical spectra from the SNLSE model tally with the experimental ones. 

It is interesting to note that the curves of Fig. la look more noisy than those of Fig. lb 
although we have performed the same number of stochastic realizations in both cases. This 
is because the number of realizations needed to achieve a given precision on the expectation 
value of a quantum operator is a function of the relative value of the spatial step h and the 
typical distance over which the nonlinear effects act; both are different in the simulations 
of Fig. and Fig. EJd. In practice, the smaller the spatial step, the fewer the number 
of realizations needed to achieve a given precision on expectation values. In Fig. ^ 50 
realizations are enough to estimate Se with an accuracy of about 1.5 dB. We also point out 
that the noise level visible at non phase-matched frequencies has no physical meaning. It 
can be lowered by averaging over a higher number of realizations. However, the number of 
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FIG. 1: MI spectra for various birefringence regime. Simulation parameters are listed in Tab. [IJ 

(a) Spontaneous S-MI in perfectly isotropic fiber in the anomalous dispersion regime. Black, dark 
gray and light gray curves correspond to a propagation length L = 500, 1000, 1500 m, respectively. 
The inset exhibits the pump spectral broadening due to the self-phase- modulation (SPM) effect. 

(b) Spontaneous V-MI in a slightly birefringent fiber in the normal dispersion regime. The pump is 
polarized along the slow axis, Stokes and anti-Stokes photons appear on the orthogonal axis. Black, 
dark gray and light gray curves correspond to a propagation length L = 16, 24, 32 m, respectively. 

(c) Spontaneous V-MI in a strongly birefringent fiber in the normal dispersion regime. The pump 
polarization axis makes an angle of 45 degrees with the slow axis. Stokes (anti-Stokes) photons 
appear on the slow (fast) axis. Black, dark gray and light gray curves correspond to a propagation 
length L = 10, 20, 30 m, respectively. 
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realizations needed to achieve an accurate estimation of Se in the non phase-matched part 
of the spectrum is usually very high. When tractable, a linear perturbation analysis will be 
less time-consuming. 

-20 




-0.5 0.0 0.5 1.0 1.5 2.0 2.5 
Phase mismatch parameter log 10 ( Ap [nr 1 ] ) 

FIG. 2: The figure illustrates the effect of an increasing birefringence on the weak-birefringence 
V-MI phenomenology. The parameters of the simulations are: Ao = 1550 nm, (3 2 = 60 ps 2 km -1 , 
7 = 2 W _1 km , Tfwhm = 100 ps, Pq = 400 W, and L = 40 m. The pump wave is polarized 
along the slow axis. The fiber beat-length was varied from 10 m to 5 cm. The values of A/3o and 
A/?i where deduced from Eq. Q14[) . The plot shows the maximum value of Se in the side bands 
as a function of A(3q. Black squares (gray squares) correspond to a single (double) peak structure 
(see the text). 



Because our algorithm permits to investigate intermediate birefringence, we have also 
studied the effect of group velocity mismatch on the transition from low to high birefringence 
limits. To our knowledge, this transition has never been fully investigated before. Fig. El 
shows the results of simulations with the same parameters as in Fig. except that the 
propagation length has been set to L = 40 m and that the fiber beat length Lb was varied 
from 10 m to 5 cm. The value of the maximal energy spectral density Se in the side bands 
is plotted versus the phase mismatch parameter. By lowering Lb, we increase the value of 
the phase mismatch A/3q and the group-velocity mismatch Aj3\ according to the relations 

Ao 



A(3 = Aft = — . 

L B cL b 



(14) 



where c is the vacuum speed of light. When the birefringence increases, the side bands move 
away from the pump spectrum and their amplitude decreases. Subsequently the side bands 
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acquire a double peak structure (see Fig. EJ). This behavior is due to the walk-off of the 
produced Stokes and anti-Stokes photons. One easily shows that l/v ga — l/v gs — a/S/^AA), 
where v gs and v ga are the Stokes and anti-Stokes photons group velocities, respectively. 
Applying this formula to our simulation and taking A/3o = 10 m" 1 , one sees that the 
Stokes and anti-Stokes photons have walked 87.6 ps away while their FWHM duration is 
100 ps. Stokes and anti-Stokes walk-off limits the coherent exponential amplification of 
quantum noise. The typical length scale over which the side bands growth takes place is 
given by T FWH m/ V^^^Po- We point out that this analysis also hold for a cw-pump: The 
coherent amplification of the side bands stops when the walk-off of Stokes and anti-Stokes 
photons exceeds the coherence length of the pump. However, in the cw-case, Stokes and 
anti-Stokes photons generated in the first coherence length act as an input noise that will 
be amplified in the following coherence length. The process is reproduced as many times 
as the number of coherence lengths in the propagation distance. One usually argues jfj, I^J 
that the weak birefringence phenomenology disappears because the coherent-coupling terms 
in Eqs. (jlH|) (those containing the factor exp[±(2)zA/?o£]) average to zero when A/3q is high. 
This statement is equivalent to saying that f2 max tends to infinity. Our analysis shows 
however that walk-off has an even stronger effect. 

Until now we have not yet demonstrated that modeling vacuum-fluctuations induced MI 
using SNLSE predicts the correct values of Se- In subsect ions II V B I we compare the absolute 
values of the energy spectral density at the maximum gain obtained using our program and 
the linear perturbation analysis. Having clarified in this way the interpretation of the results 
of our numerical simulations we turn to a detailed comparison of the effect of classical and 
quantum noises. That is we compare the effects of classical noise in the initial conditions 
and the quantum noise added at each step of the integration. 

B. Comparing numerical integration and linear perturbation theory 

In our numerical simulations we have taken the pump laser to be a Gaussian pulse without 
chirp. Its instantaneous power and energy spectral density can be written 



Pit) 



P exp(- — 2), 



(15a) 



t 




(15b) 
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with a t cr^ = \- The numerical integration of the SNLSE provides us with the spectral 
density of energy Se{L,Q) at the end of the fiber, see Eq. (|12p. On the other hand the 
linear perturbation theory is based on small perturbation analysis around a continuous 
monochromatic pump. We would like to compare quantitatively the predictions of these 
two approaches. 

For definiteness we carry out this comparison in the case of scalar MI. We shall focus our 
investigation on the intensity of the sidebands at the peak of the MI gain (Q = f2 max ) in the 
two approaches when we modify the propagation length and as we modify the duration <j t 
of the pulse. 

The linear perturbation theory is based on a continuous monochromatic pump. For this 
reason the theory predicts a rate of photon production per unit time. This suggests that if 
one takes the pump to be a pulse localized in time, the number of photons produced should 
be proportional to the pulse duration, all other parameters being kept constant. Simulations 
based on SNLSE confirm this phenomenology. This is illustrated in Fig.Ekj where the energy 
spectral density at f2 max is plotted as a function of ^PqL for three different pulse durations. 
Note however that the above argument is valid for square pulses but is not very satisfactory 
for Gaussian ones. A better understanding of the origin of this scaling can be obtained by 
making appeal to the notion of mode and of Heisenberg box. 

In order to introduce the notion of mode, recall that a temporal signal can be repre- 
sented by a distribution in the time-frequency plane. But because of the time-frequency 
uncertainty relations, a point in this plane has no physical meaning. This problem is well 
known in signal processing where one usually thinks in terms of local time-fre quen cy de- 
compositions, using windowed Fourier transforms (WFT) or wavelet transforms 2M- Such 
a local time-frequency decomposition allows one to decompose a signal into orthogonal local 
functions, called modes. These can be depicted as surface elements in the time-frequency 
plane. Fourier-transform limited pulses, such as our pump pulse (Jl5|) . are represented by 
Gaussian distributions on the time-frequency plane, called the Wigner-Ville distributions. 
This distribution is very similar to the Wigner distribution used in quantum optics to rep- 
resent a quantum state of light. In particular two different Fourier-limited Gaussian pulses 
sufficiently different in time or central frequency can be considered as quasi-orthogonal 
modes. A set of quasi non overlapping Gaussian Wigner-Ville distributions can be taken 
as a base for the time-frequency decomposition of the field. One can visualize this modal 
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FIG. 3: Values of the number of photons created at the maximum gain frequency obtained by 
integrating numerically the SNLSE as a function of propagation length L. In order to keep the 
maximum gain frequency constant we have kept the peak power Pq constant in each figure. The 
horizontal axis is given in dimensionless units of jPoL where Pq is either the power of the continuous 
pump wave in the linear approximation, or the peak power of the Gaussian pump pulse. The top 
panel is plotted in the density of energy representation whereas the bottom panel is plotted in the 
number of photons-per-mode representation using the rescaling of Eq. Q16JI . In both panels, the 
up-triangles, the circles and the down-triangles correspond to FWHM durations respectively equal 
to 4 ns, 1 ns and 0.25 ns. Note that these three curves coincide perfectly in panel (b), thereby 
showing the relevance of the rescaling H16j) . In the bottom panel, the dash-dotted line results from 
the analytic solution Eq. @. The dashed line corresponds to Eq. Q convoluted with the pump 
shape then rescaled according to Eq. ((TB|l . see the text. All other parameters are identical to those 
used in Fig. 
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base by imagining that the time- frequency plane is paved with adjacent elementary surface 
elements, called Heisenberg boxes, roughly representing the area of Gaussian quasi non over- 
lapping Wigner-Ville distributions. The precise area of the Heisenberg boxes is a matter of 
taste depending on how strong orthogonality is required. A usual convention is to take this 
area equal to a t x a w = 1/2. 

This set of (quasi) mode is convenient for our problem because, during the pulse propaga- 
tion in the fiber, the uncertainty on the creation time of a photon is defined by the variance 
al of the pump pulse, and implies an uncertainty on its frequency defined by the variance 
cr^j. More precisely, in the case of S-MI, the pump pulse (fTHj) produces photons that occupy 
single Heisenberg boxes located at the same time as the pump but at angular frequencies 
u ±Q (0 < Q < v^^max)- Formulas (jHJ) and (jHJ) of Sec. Ill Bl thus give the number of photons 
created in these Stokes and anti-Stokes modes. 

Our numerical simulations provide us with the spectral distribution of energy Sg. Hence 
we need to reexpress this as the number n of photons produced per mode of duration a t and 
spectral width <j w : 

n(n) = xa LU = -r^- • (16) 

niOQ moo i fwmh 

Fig. Eb shows the effect of scaling the spectra according to Eq. (fTfijl. When expressed in 
terms of number of photons per modes the three curves of Fig. EH (corresponding to three 
different pump durations) come down to a single one in Fig. Eh (continuous line). This shows 
that the notion of mode helps in interpreting the results of the numerical integration. The 
number of produced photons does not depend of the pump duration: The pump duration 
just alters their time-frequency characteristics. The dash-dotted curve in Fig.Eb corresponds 
to a direct application of Eq. (JHJ). A discrepancy with the simulations based on the SNLSE 
can be noted. It is simply due to the fact that our choice of size of Heisenberg boxes, hence 
the normalization factor in Eq. (|16|h is somewhat arbitrary. By taking the Heisenberg boxes 
a bit bigger, one can put the continuous and dash-dotted curves of Fig Eh in superposition. 
In the rest of the text, however, we maintain the normalization relation Eq. (|T6|) for clarity. 

There is another way to compare the results of SNLSE simulations to the linear per- 
turbation theory that avoids the concept of modes. One first computes the power spectral 
density given by the quantum perturbation analysis on a monochromatic pump wave. Then 
one convolutes this with the Gaussian spectral distribution of the real pump pulse. This 
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procedure gives a good approximation of the energy spectral density generated by the Gaus- 
sian pump. The dashed curve of Fig. EJd corresponds to the peak energy spectral density 
computed by this method and rescaled according to Eq. in order to be independent 
of the pulse duration. The agreement with the continuous curve is now much better. The 
origin of the small difference in slope in the exponential amplification regime still remains 
unclear. It may be due to the self-phase-modulation broadening of the pump spectrum, 
which is not taken into account by the linear perturbation analysis. 

As a conclusion we obtain a good quantitative agreement between the simulation based 
on SNLSE and the linear perturbation analysis. We have also shown that the only physical 
quantity that can be rigorously predicted by the quantum nonlinear propagation theory is 
the energy spectral density and that the concept of mode, although very useful, must be 
handled with care. Especially formulas like (JBJ) and/or Q may only be used as a rough 
approximation tool because there is no objective way to define a time-frequency mode. 

C. Classical versus Quantum Noise 

From Figs. H3 it is clear that the spontaneous MI growth can be divided into two different 
stages. So long as the number of particles created by mode is less than one (n < 1), one 
is in a quantum regime dominated by vacuum fluctuations. In contrast, for n above 1 the 
modulation instability is amplified exponentially and quantum effects become negligible. 

We now investigate the transition between the quantum and classical regimes in presence 
of some classical noise. We have chosen to model this noise by modifying the initial condi- 
tions and adding a white noise to the amplitude of the pump pulse in the Fourier domain. 
For definiteness and simplicity, we continue to focus on scalar modulation instability. 

The classical initial noise N(Q) was chosen according to two criteria: (i) the noise must 
correspond to a random fluctuation of the pump amplitude in the time domain and (ii) the 
statistic of N(Q) (for each frequency) must lead to a spectral density of energy (N(Q)N*(Q)) 
constant as a function of the frequency. Several noise definitions can meet these criteria. 
We have chosen to study two particular cases : the pure spectral phase noise 




(17) 
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and the Gaussian noise _ 

N G (n) = A x ( Cl (n) + i C2 (n)), (is) 

where A is a real constant and Cj{ty are independent real zero-mean Gaussian white noise 
random fields. In our simulations £j(^) where discretized and replaced with random quan- 
tities (one for each discretized frequency) drawn according a zero-mean Gaussian law 
of variance 1. We have compared S-MI spectra averaged over 50 realizations, obtained from 
both classical noises (with quantum noise set to zero). Both noises lead to equivalent results. 
The difference between the spectra is lower than 0.3 dB on the full spectral span, which is 
less than the residual averaging noise (See FigJIJi), and therefore negligible. Hereafter, the 
classical noise is set according to Eq. (|17|) . It corresponds to a white flux of photons without 
any phase correlations between each frequency component. 

In Fig. |3] we compare the peak intensity of side bands in the case of vacuum- fluctuations 
induced MI and the (unphysical) case of MI induced by the classical noise alone. In the 
quantum regime where the number of particles per mode is small the two approaches differ 
strongly whereas in the exponential amplification regime they give similar results, although 
different classical noise levels give rise to different final number of photons. 

The simulations of Fig. take into account both classical and quantum noises. They 
illustrate the realistic situation when both noises compete for producing modulation insta- 
bility. If the order of magnitude of the classical noise is such that there is less than one 
photon per mode, the quantum noise dominates and the curve is (except for very small 
values of the gain) indistinguishable from the case where there is no classical noise. On the 
other hand if the number of noise photons per mode is much larger than one, the classical 
noise dominates and the intensity of the sidebands is indistinguishable from purely classical 
situations depicted in Fig. HJ These numerical results are consistent with Eqs. (jSJ), in which 
the term 1 represents the contribution of the vacuum fluctuation. If n s or n a are higher than 
1, Eqs. ((HI) tend to Eqs. flU) corresponding to the classical description. 

In summary we have shown that if one considers only the peak intensity of the sidebands, 
the quantum origin of the instability can only be seen in the regime where the number of 
photons per mode (produced or initially present) is small whereas when the number of 
photons per mode is large the effect of vacuum fluctuations is indistinguishable from that of 
classical noise. Good quantitative agreement between the two approaches in the exponential 
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FIG. 4: Comparison between a MI growing from quantum noise (solid line) and from purely classical 
noise. The number of photons in the Stokes and anti-Stokes modes is plotted as a function of the 
propagation length (expressed in dimensionless units ^PqL). The input noises in the classical 
situation correspond to the following amount of photons per mode: (i) 1/40 (dash-dotted line), (ii) 
1/2 (dotted line), and (iii) 10 (dashed line). The simulation parameters are those of Fig. 

amplification regime is obtained when the number of noise photons is 1/2 per mode. Note 
that other quantum effects, such as two mode squeezing, may be present in the regime where 
many photons are produced per mode, but exhibiting them requires looking at correlations 
between the two sidebands. 

D. Using classical noise to compute the instabilities induced by vacuum fluctua- 
tions 

We can now discuss the well-known trick which consists in introducing a half photon per 
mode into a classical simulation (then removing it) to simulate the spontaneous effects. To 
see how this works we compare Eqs. (J3J which describe the MI induced by classical noise 
and Eq. (jHJ) which is derived from the quantum theory and gives the number of photons 
produced per mode by the action of vacuum fluctuations alone. Now, if we introduce the 
same amount of Stokes and anti-Stokes noise photons no in Eqs. (0J) we find: 



n^ 1 1 rig 1 



1 



sinh 2 (7P £) 




(19) 
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FIG. 5: Simulated values of the number of photons in the Stokes and anti-Stokes modes in the 
realistic case when MI grows from both quantum and classical noise. The classical noise intensities 
are defined in order to correspond to an amount of photons per GHz: (i) 0.1 (circles), (ii) 1 
(up-triangles), (iii) 10 (down-triangles), and (iv) 100 (stars). The squares correspond to a purely 
quantum noise. The simulation parameters are those of Fig. ^i. The inset is a zoom corresponding 
to the the lowest 7-Po-^-P arame t ers - It is drawn in a logarithmic-scale. 

where CI and Qu denote respectively the classical approach and the quantum approach. 
Taking uq — h, Eq. (|T§|) shows the agreement between the quantum predictions (right hand 
side) and the classical predictions (left hand side). Note that n can be any real value 
(except for 0), hence the spontaneous growth of the number of photons per mode can be 
simulated with any number of initial classical photons if the pump depletion is neglected. 

We have compared, using numerical integration, the direct quantum approach based on 
the SNLSE (without classical noise) and the classical approach in which one first integrates 
the NLSE with some initial classical noise then rescales the spectra according to the left side 
of Eq. (fT£i|). The discrepancy between both approaches is measured by the following ratio 
in dB scale: 



Applying Eq. (f2*U|) to the data reported in Fig. 01 one founds that for ^PqL higher than 0.2, 
7] is constant for any classical noise amplitude: 77=1.9 dB ±0.3 dB. Note that both noise 
definitions Eq. (fTTj) and Eq. (fTHj) lead to the same rj. Moreover this results may be extended 
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to ^PqL lower than 0.2 by increasing the number of realizations. The origin of non-zero value 
of rj may come from the mode definition used in Eq. (|16p. Indeed the numerical integration 
of SNLSE gives the physical spectral density of energy Se whereas the calculus trick — left 
hand side of Eq. (|19|) — leads to values interpreted as a number of photons per mode which 
must be rescaled to give Se- 

In summary, in the case of pump pulses of finite duration, a full quantum treatment based 
on the SNLSE leads to a direct quantitative prediction whereas the calculus trick only gives 
a good approximation whose accuracy is dependent of the modes definition. 



V. CONCLUSION 



We generalized to an arbitrary level of birefringence the stochastic nonlinear Schrodinger 
equations describing the propagation of pulses through a nonlinear medium with linear 
birefringence and group-velocity dispersion, and developed numerical routines to compute 
them. Because these stochastic equations are equivalent to quantum field operator equa- 
tions, we used them to compute spontaneous (or vacuum-fluctuations induced) modulation 
instability spectra in various birefringence regimes, including weak, high but also intermedi- 
ate birefringence which has not been studied so far. In particular we showed that the decline 
of the number of photons generated by the weak-birefringence V-MI when the birefringence 
increases, is attributable to the increase of the walk-off between Stokes and anti-Stokes pho- 
tons, although the weak-birefringence V-MI gain remains constant. We then investigated 
the absolute values of the energy spectral density at the maximum gain in the case of scalar 
modulation instabilities induced by vacuum fluctuations. We obtained good quantitative 
agreement between the simulation based on SNLSE and the linear perturbation analysis. 
Then we have carried out a detailed comparison of the effect of classical and quantum noise 
and shown that the quantum origin of the instability can only be seen in the regime where 
the number of photons per mode produced or initially present is small. Finally we note that 
the quantum nonlinear propagation theory predicts the energy spectral density and that the 
concept of mode, although very useful, must be handled with care. The present work forms 
the basis for numerical and experimental investigation of vacuum-fluctuations induced V-MI 
in regimes which have been little investigated so far, and we hope to report on this in the 
near future 131. 
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Although we have not developed this aspect in this article the stochastic equations can 
also be used to computed intensity correlations between side bands and predict special 
quantum effects like squeezing. For this reason the stochastic equations (|13|) are a valuable 
tool for computing quantum effects in birefringent nonlinear media, especially optical 
fibers. The stochastic model may also be adapted to include Raman and Brillouin effects 



sec 



12j for the scalar case). Higher order dispersion effects can also be included in a 



straightforward way. 
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APPENDIX: STOCHASTIC COUPLED NONLINEAR SCHRODINGER 
EQUATIONS 



In order to derive the SNLSE (|13jh we will proceed in three stages. First, we will estab- 
lish the interaction Hamiltonian that governs field evolution in a lossless, dispersive, and 
birefringent fiber. We will only present an heuristic derivation of this Hamiltonian and put 
the emphasis on appropriate approximations. A rigorous derivation requires a discussion of 



electromagnetic field quantization in material media [12J, 25, 26 1, which is outside the 
scope of this article. Second, we will use this Hamiltonian to find the Liouville equation de- 
scribing the evolution of the density operator of the field and convert it into a Fokker-Planck 
equation using the p( + )-representation. Finally, we will establish the connection between 
the Fokker-Planck equation and the stochastic equations (fT^j) . 
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1. Linear and Nonlinear Hamiltonians 



In a dispersive birefringent medium the positive-frequency part of the s-polarized electric 
field component (s = x, y) can be decomposed on monochromatic modes in the following 
way 0: 

^-iJ^^f^F^. (A,) 

Eq. ()A.1|) can be seen as defining a s {(3). The operators a s {[3) and its hermitic conjugated 
a\{(3) are, respectively, the annihilation and creation operators for a s-polarized photon 
propagating in the fiber with a propagation constant (3 and having an angular frequency 
u> s (/3). According to the canonical quantization, they satisfy the commutation rule 

[a(P),tf{J?)] = 5(P-p). (A.2) 

In Eq. flA.lj) . n s (j3) and v gs ((3) are respectively the linear index of refraction and group 
velocity corresponding to the s-polarized monochromatic mode with frequency u) s (/3), and 

A= \F(x,y)\ 2 dxdy. (A.3) 



The operator representing the s-polarized electric field component is given by 

E s (r) = Ei + \r) + Ei-\r), (A A) 

where Ei \r) = [Eg(r)]^ is the negative-frequency part of E s (r). 

The total Hamiltonian representing the sum of the vacuum electromagnetic energy 
and the dielectric energy stored in the fiber can be decomposed in a linear part Hi and a 
nonlinear one H^i. 

The linear part, 

H L = fdpYl t^s{P)a\{P)a s {(3), (A.5) 

s=x,y 

takes into account the free-field energy and the energy stored into the dielectric through 
linear interactions, including the effects of linear dispersion and linear birefringence through 
the dispersion relations uj s = ui s {f3). 

If the field bandwidth is narrow compared to the central angular frequency Uq, dispersion 
can be neglected in the interactions and the nonlinear part of the Hamiltonian can be 
written 

H NL = -\e d J d 3 r XmE^E^E^E^, (A.6) 

ijkl 
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where Xijki stands for Xijkii^o] uq, — uq). This simplified Hamiltonian only takes into 
account the Kerr effect, which is the dominant one for quasi- monochromatic fields. Since 
the medium is supposed lossless, the x^ tensor has the full permutation symmetry 
The degeneracy factor d = 6 takes this symmetry into account, by counting the number of 
way to permute the frequency arguments and the indexes of the x^ tensor. A further useful 
approximation is to consider that the x^ process is isotropic |27l . l28j | : 

Xijkl Xxxyy^ij^kl Xxyxy^ik^jl Xxyyx^il^jk- 

Since the field bandwidth is supposed narrow, the permutation symmetry also requires that 

Xxxyy = Xxyxy 

"Xijkl Xxyxyi&ij&kl &ik&jl) Xxyyx&il&jk- i.-^-^) 

Using Eq. (jA.7j) . the nonlinear Hamiltonian becomes 

Hnl = -\^Xxxxx I d\ ( ^E^E^E^E^ + (1 _ a) ££<->£H£M£« 

+BJ2E^E^E^E^), (A.8) 

s^s' ' 

where we defined B = Xxyyx/Xxxxx, and factored out Xxxxx = ^Xxyxy + Xxyyx- 

Another consequence of the narrow-bandwidth assumption is that the square-rooted 
bracket in Eq. (jA.l|) can be taken out of the integral and one can write 

£(+) (r) w i f ^' | 1/2 f(^)i( Z| t)e'^, (A.9) 
V2e n s0 c A) 

where 

i(uj t-f3 s0 z) r 

i> s (z,t) = -=— / d(3 s a s (/3 s ) e^ z . (A.10) 

V 2vr J 

In Eqs. (lA~9l) and (lA~IOl) . 

n so, v gso, and [3 s q stand respectively for the index of refraction, the 
group-velocity, and the propagation constant at frequency u on the s-axis. The operator ip 
is an envelope operator because fast oscillations in space and time have factored out. This 
implies that ?/> is explicitly time-dependent in the Schrodinger picture. The operator iplip s dz 
represents the number of s-polarized photons in [z, z + dz\. One can easily check that 

$ a (z,t),$ a ,y,t)] = 5 aa ,8(z-2>). (A.ll) 

Using Eq. (jA.9|) . H NL takes the following simple form: 
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Hnl = ~Q 



ft x ftJJ x + ^li>yi>y) + 2(1 - B)^Jl^ y 

+B((&^y iA ^ + (^ 2 x e-^)]dz, (A.12) 
where 



3hwy goXxxxx a 2 

Let's note that we have set uq = n X Q « n^o and t> 9 o = f ga ;0 ~ ^o- The linear Hamiltonian 
Hl can also be expressed in function of the operators (ip s , fy), s = x,y, by developing u) s {(3) 
in a Taylor expansion around (3 s q up to the second order, 

u 8 (J3) =lu + ^(/3 - (3 s0 ) + y^~ ^ + ^ A - 14 ) 

where u' s = %|/3 = f 9S o and cu" = ^flftr Using Eq. ()A.14|) and inverting Eq. ()A.10|) . one 
finds that Hl = U + where 

£> = hw J $(z)i> s (z)dz, (A.15) 

s=x,y 

and 

H L=2 22J + w --d7-d7 ]daf - (A - 16) 

In the Heisenberg picture, the hamiltonian C7 is responsible of a free oscillation exp(— ic^ot) 
of the fields ^ (s = This oscillation will cancel out the explicit oscillation exp(iu;ot) 

already present in the definition (|A.10J) . For this reason, we prefer to continue the discussion 
in the interaction picture: 

(ips)i = i> s exp(-iuJot), (A. 17) 

Hj = H T -U = H' L + H NL . (A. 18) 

To simplify the notations we will drop the I index in later equations. 
2. From Liouville to stochastic equations 

In the quantized theory, the state of the electromagnetic field is represented by the density 
operator p(t). Its evolution, in the interaction picture, is given by the Liouville equation 

ih^p=[H,p], (A. 19) 
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where H is the Hamiltonian defined by Eq. (jA.18|) . Using Eqs. (jA.16|) and ()A.12|) . the 
calculation of the right-hand side of Eq. (jA.19|) is straightforward, so we do not write it here 
explicitly. 

In order to obtain stochastic equations from the Liouville equation (jA.19|) . we generalized 



221 for monomode fields and their extension to 



the argument of Drummond and Gardiner 

mu , lm „de sea,, Se.ds gl ve„ is Q. We ia,o d uee the m ulti mod e eoherent states ,{„}> 
defined as the eigenstates of the annihilation operators d s (/3) 

a s (/?)|{a})=a s (/?)|{a}>. 

As a consequence, \{a}) are also eigenstates of the envelope operator (|A.17j) 

^ s (z)\{a})=^ s (z)\{a}), 

with 

i/>.(z) = -L ! d(3 s a s ((3 s ) e *03.-A°)*. 



/2tt 

This suggest the alternative notation \i^(z)) = |{a}), with i^(z) = (i[) x (z), if) y (z)), for the 
multimode coherent state. The basic idea of P( + )-representation is to expand the density 
operator on nondiagonal coherent state projection operators defined as 

M*(z)) \^M^m\ (A 20) 

where ip\z) = (ip].(z),^(z)) is a new set of fields different from if){z). Denoting &(z) = 
(ip x (z),ipl(z),ipy(z),il>^(z)), this expansion can be written in the following way: 

p(t) = J P(*;t)A(*)d//(tf), (A.21) 

where the integration measure dfj,(^) means that the integration is carried over all the 
possible fields ijj s and s = x, y. Taking into account the definition (|A.20|) one can show 
that, 

40)A = i> s {z)k (A.22a) 
ft(z)k = L\(z) + T1 ^) A, (A.22b) 



6ip s (z) 

Ut(z) = ^(z)k, (A.22c) 

= (tTTTT + M*)) A > ( A -22d) 
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where 8/8ip s (z) and 5/5ipl(z) are functional derivatives. Eqs. (| A. 22(1 generalize the corre- 



sponding monomode identities of 
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The P-function always exist and is positive for any density operator. The P-function is 
useful for calculating normal ordered moments: 



<($) m (Vv) n > = / M m fo) n P(*',t)dn(?) 



(A.23) 



In particular, Eq. ()A.23|) shows that P is normalized to unity: 1 = J P(^;t)dfi(^). It can 
be interpreted as a genuine probability density on the (infinite dimensional) space sustained 
by the field ip s (z) and ipl(z). 

To obtain the time evolution of P, we insert the expansion (jA.21|) into the Liouville 
equation (jA.19|) and find 



?^A(#)dp(*) 



AW 



(A.24) 

where summation over k and / is implied. The CVs are the components of a four- dimension 
drift vector C, with 



— ^ 



+ z- 



x dz 2 

+ieB4>itp 2 y e- 2iAf3oz 



(A.25) 



The C*2, C3, and C4 components have a similar form. C2 is obtained from (jA.25|) by making 
the substitution (i) % — > — z, ^ <-> and if> y i?!. C 3 is obtained by (ii) exchanging x- 



and ^/-indexes in (1A.25J) . and making the subtitution A/^o — >■ — A/3q. To obtain C4, both 
substitutions (i) and (ii) must be performed. The Dki are the elements of a symmetric 
diffusion matrix D that can be written in the form D = BB T ', where 



/V* 

i^t 
% 



B^ y e- tAI3oz 



\ 

iy/Bi^e^* 
B^ x e lA,3oz 

-i\fB^\.eT iA ^ z ) 



(A.26) 



Using Eqs. (lA~22l . on can deduce from Eq. (jA.24|) that the P(\&; £) verifies a functional 

n 

Fokker-Planck equation with a semi positive-definite diffusion matrix. We refer to [9] for a 
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demonstration since Eq. (jA.24|) has the same structure as Eq. (4.19) in [9|. Because of the 
semipositivity of the diffusion matrix, the positivity of P is maintained during evolution. 

The stochastic equations equivalent to the Fokker-Planck equation for P can be written 
in the following compact form: 



where (k, I) e {1,2,3,4} 2 , and Q(z,t) are the independent zero-mean Gaussian white noise 
random fields introduced in Sec. IIHI and characterized by the second order moments (fTT| . 
The C vector gives the deterministic evolution of the fields as predicted by the classical 
theory of light. The way vacuum fluctuations modify the classical evolution is determined 
by the structure of the B matrix. If one discards the stochastic terms, the fields ip\ and ip s 
appear to be just complex conjugated of each other. However, when vacuum fluctuations are 
taken into account, ip\ and ip s must be treated as independent fields that are only complex 
conjugate in mean. 

As they stand, Eqs. (jA.27|) seems to differ from Eqs. (fT3"j) . Actually, both writings are 
equivalent. To highlight the equivalence we first notice that, according to the instantaneous- 
power normalization of the (A s , A\) fields, one has the following relations: 



d_ 

Of 



y k (z,t)=C k (*) + B kl (V)( l (z,t), 



(A.27) 



A s (z,t) = y/fiu} Vg a0 1p 3 (z,t), 



(A.28a) 



A\(z,t) = ^hw v gs0 ipl(z,t). 



(A.28b) 



Inserting (jA.28|) into (jA.27|) . and noting that uj' s = v gs0 and 6 = hu} Vg 'y, we find 
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dA x ^ 1 dA a 



dz v gx0 dt 



dAl ^ 1 dAl 



+i ^flt X + *7 + (1 " B)A\A y ] A x + z 7 5(^) 2 4e- 2 ^ 2 

(A.29a) 



'2v gx0 dz 



< d 2 Al 



dz v gx0 dt 



dA y ^ 1 dA y 

dZ VgyQ dt 



dAl l dA], 
+ 



—% 



2v gx0 dz 2 
u" d 2 A 



+^ 



y " ^y 



dZ VgyQ dt 



— I 



- i 1 [A X A X + (1 - B)A\A y ] 4 - z 7 5(4) 2 ^ e+2lA/3 ° 2 

(A.29b) 

+ i 1 [A\A y + (1 - B)A\A X ] A y + i 1 B{A x ) 2 Ale +2lA ^ z 

(A.29c) 

- %1 [A\A y + (1 - B)A\A X ] A\ - i^BiA^Aye- 2 ^ 

(A.29d) 



2v gy0 dz 2 



2u 9y0 <9z 2 



Eqs. (| A. 29(1 differs from Eqs. ()13|) only in the first term of the right member. For each axis, 
the group- velocity dispersion parameter is /3 2s = — u"/Vg s0 . When the typical pulse duration 
T is such that T/fas is much bigger than the group- velocity v gs o, which is the common 
situation in fiber-optics, the following operator approximation holds 



d 2 



l d 2 



dz 2 v 2 s0 dt 2 ' 



(A.30) 



Inserting (|A.30J) into Eqs. (|A.29J) . and noting that usually (3i x ~ foy = one obtains the 
stochastic equations (fT3Jl . 
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